Exotic phase diagram of a cluster charging model of bosons on the kagome lattice 
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We study a model of hard-core bosons on the kagome lattice with short-range hopping (t) and 
repulsive interactions (V). This model directly maps on to an easy-axis S = 1/2 XXZ model on 
the kagome lattice and is also related, at large V/t, to a quantum dimer model on the triangular 
lattice. Using quantum Monte Carlo (QMC) numerics, we map out the phase diagram of this model 
at half-filling. At T = 0, we show that this model exhibits a superfluid phase at small V/t and an 
insulating phase at large V/t, separated by a continuous quantum phase transition at V c /t « 19.8. 
The insulating phase at T — appears to have no conventional broken symmetries, and is thus a 
uniform Mott insulator (a 'spin liquid' in magnetic language). We characterize this insulating phase 
as a uniform Z2 fractionalized insulator from the topological order in the ground state and estimate 
its vison gap. Consistent with this identification, there is no apparent thermal phase transition upon 
heating the insulator. The insulating phase instead smoothly crosses over into the high temperature 
paramagnet via an intermediate cooperative paramagnetic regime. We also study the superfluid-to- 
normal thermal transition for V <V C . We find that this is a Kosterlitz-Thouless transition at small 
V/t but changes to a first order transition for V closer to V c . We argue that this first order thermal 
transition is consistent with the presence of a nearby Z2 insulating ground state obtained from the 
superfluid ground state by condensing double vortices. 

PACS numbers: 75.10.Jm, 05.30.Jp, 71.27. +a, 75.40.Mg 



I. INTRODUCTION 



Spin liquids are quantum disordered paramagnetic 
phases that preserve all lattice symmetries. While there 
has been considerable progress in understanding the ef- 
fective field theories and properties of such spin liq- 
uid phasesi^^^, showing that the excitations in this 
phase carry fractional quantum numbers and interact 
with emergent gauge fields, there are very few micro- 
scopic models which can be convincingly shown to ex- 
hibit a spin liquid phase and they fall, roughly, into two 
categories. One class of microscopic Hamiltonians are 
quantum dimer models^, which have been proposed to 
describe spin gapped phases of quantum magnets. Some 
of these models, on the triangular and kagome lattice in 
two dimensions, exhibit Z2 fractionalized quantum disor- 
dered phases. However, the Hilbert space of such quan- 
tum dimer models has a strong local constraint, namely, 
the number of dimers emerging from a site is fixed. It is 
therefore interesting to examine other models which do 
not have such local Hilbert space constraints and, thus, 
no extra conservation laws other than the total spin or 
total S z . Under this second category are models, which 
incorporate a so-called "cluster charging" energy^, which 
penalizes spin configurations where S z summed over a 
local "cluster" of sites differs significantly from its mean 
value. Some simple models of this type can be shown to 
reduce to effective quantum dimer Hamiltonians in the 
limit of a large charging energy, with the Hilbert space 
constraints emerging at low energy from energetic consid- 
erations. By the usual mapping between 5=1/2 spins 
and hard core bosons, such quantum spin models can 
be alternatively viewed as boson models. The spin liq- 
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FIG. 1: (color online). The schematic phase diagram of the 
model given by Eq.(l) or equivalently Eq.(2). The superfluid- 
insulator transition in the model (1) corresponds to the XY 
ferromagnet to spin liquid transition in the model (2). No- 
tice that the Z2 fractionalized Mott insulator (Z2 spin liquid) 
exists only at zero temperature in two dimensions and as a 
consequence, there is no finite temperature phase transition 
in the paramagnetic region, but only the crossover between 
different regimes as explained in the text. Near the quantum 
critical point between the superfluid and the Z2 fractionalized 
insulator, there should be a quantum critical region (denoted 
as QC in the figure) that is not discussed in the present paper. 



uid phase of the spin model at zero magnetization then 
corresponds to a uniform Mott insulator of bosons at 
half-filling. 

In this paper we study, using a generalized stochastic 
series expansion QMC algorith m 10 ! 11 , a "cluster charg- 
ing" model of hard core bosons on the kagome lattice 
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with the Hamiltonian 

Here b\(bj) is the boson creation(annihilation) operator, 
t > is the hopping amplitude, V > is the repulsion 
strength, = is the number operator, and /i = 12V 
is the chemical potential which fixes the boson density 
to be at half-filling. The hopping term connects only the 
first, second and third neighbors. The repulsive interac- 
tion is similarly short-ranged. The main result of this 
paper is the phase diagram summarized in Fig. ([Tj). 

We begin by reviewing our earlier results 1 ^ which show 
that model (|TJ) exhibits, at T = 0, a superfluid-insulator 
quantum phase transition at (V/t) c « 19.8. The insulat- 
ing phase is a uniform, topologically ordered, Mott 
insulator. We then present new results for the finite 
temperature phase diagram. We find that the insula- 
tor crosses over into the high temperature normal phase 
via a cooperative fluctuation regime but with apparently 
no intervening thermal phase transitions. This bolsters 
the case for no broken lattice symmetries in the insulat- 
ing phase. On the superfluid side, for small values of 
V/t, raising temperature leads to a Kosterlitz-Thouless 
(KT) transition from the superfluid phase to a normal 
phase. As we approach the quantum critical point how- 
ever, with V/t > 13, the KT transition converts into a 
first order transition. We argue, from a renormalization 
group analysis of an appropriate sine Gordon model, that 
this is consistent with an underlying quantum phase tran- 
sition into a Zi insulating phase driven by double vortex 
condensation. 

Our model in Eq. (fT]) is inspired by an XXZ spin model 
for S = 1/2 quantum spins proposed by Balents, Fisher 
and Girvii>i2. (BFG), with a Hamiltonian 

tfxxz = - Jl £[(S£) 2 + (S y Q ) 2 3] + J z J2( S o) 2 - (2) 

o o 

Here Sq = Xaeo &i 1S a sunl over * ne s ^ x s P ms on eacn 
hexagon of the kagome lattice unit cell, denotes a 
sum over all hexagons on the lattice. This model is eas- 
ily seen to be a short-range anisotropic XXZ model, with 
only the first, second and third neighbor interactions be- 
ing nonzero and equal to each other. 

Analyzing the model in Eq.© for J± < and 
Jz/\J±\ ^> 1, and interpreting each up-spin as a dimer 
on a triangular lattice, BFG showed 1 ^ that the Hamil- 
tonian is dual to an effective triangular lattice quantum 
dimer model with three dimers touching each site. In spin 
language, this effective model takes the form of a ring- 
exchange model, with an exchange scale </ r i ng = J\/J z , 
which describes quantum dynamics in the Hilbcrt space 
with Sq = on each hexagon, the local constraint 
arising from energetic considerations at large J z /J±. 
Supplementing this model with an additional four-site 
(Rokhsar-Kivelson (RK) 14 ) potential term of strength V4 



they showed that this modified Hamiltonian is in spin liq- 
uid phase for V4 = J r ing, which was argued to be stable 
for small deviations V4 < J T i ng . Later exact diagonal- 
ization (ED) numerics 1 ^ showed that the ring-exchange 
model appears to be in this spin-liquid phase down to 
t>4 = 0, but only system sizes upto 20 unit cells could be 
explored. 

The relation of the XXZ model to the hard core boson 
model we study follows upon using the standard mapping 
between 5=1/2 quantum spins and hard core bosons. 
Specifically, the Hamiltonian we study in this paper at 
half-filling is equivalent to that in Eq. © if we set J± = 
t > and J z — V > 0. Since the ring-exchange physics 
is independent of the sign of Jj_, we expect to recover, 
for large values of J z , the physics of the XXZ model with 
Jj_ < studied by BFG. On the technical side, the choice 
of J± > eliminates the QMC sign problem and allows 
us to go significantly beyond earlier work on this class of 
models. 

This paper is organized as follows. Section II reviews 
some of our earlier results and presents some new re- 
sults on the zero temperature phase diagram, including 
a discussion of topological order in the insulator. Section 

III discusses the finite temperature region of the insu- 
lating phase, including the temperature dependence of 
the energy and an estimate of the vison gap. Section III 
discusses numerical and analytical results for the finite 
temperature superfluid-normal phase transition. Section 

IV presents a summary of the results. 



II. ZERO TEMPERATURE PHASE DIAGRAM 

We begin by reviewing the zero temperature phase di- 
agram of model (jTJ) which was studied by us in an ear- 
lier paper—. For V/t = 0, the bosons only experience 
the hard-core constraint, and therefore condense into a 
superfluid phase. As we turn on interactions, the local 
charging energy V penalizes those configurations where 
the total number of bosons on any hexagon deviates from 
its mean value of no = 3. At large V/t, this cluster 
charging energy leads to locally incompressible hexagons. 
This suppresses off-diagonal long range order (and super- 
fluidity) and drives the system into an insulating phase. 
In the following subsections we review our earlier nu- 
merical results which show that the superfluid insulator 
transition is a continuous quantum phase transition with 
2 = 1. We also review our results and present new nu- 
merical data which show that the insulator has four de- 
generate, topologically distinct, ground states. 



A. Quantum superfluid-insulator phase transition 

For small values of V/t, we expect the ground state 
of model (JTJ) to be a superfluid. We characterize this 
superfluid phase by its superfluid density p s measured 
through winding number fluctuations^ W ai 2 in each of 
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FIG. 2: (color online). Upper panel: Finite size scaling of 
p s for P/L — 3/ (it). Lower panel: Data collapse of the su- 
perfluid density p s for P/L = 3/(4t), (V/t) c = 19.80(2), and 
u = 0.67(5). 

the two lattice directions, with 



(Wl) + (W ( 



2f3t 

where (3 is the inverse temperature. For small V/t, p s is 
large and its value agrees with the mean field estimate 
obtained by Gutzwiller projecting a free Bose condensate 
so as to satisfy the hard core constraint^. p s decreases 
with increasing V/t, eventually vanishing for V/t > 20 
suggesting a phase transition to an insulating phase (X*). 
This behavior is in sharp contrast to a nearly identical 
kagome lattice model where the hopping and repulsive 
interactions are restricted to the nearest neighbor only — 
in that case^ a uniform superfluid persists for arbitrarily 
large V/t. The absence of a jump in p s on going through 
the transition suggests that the ST — X* transition is 
continuous. 

In the vicinity of a continuous QPT, the superfluid 
density should scale as 



Pa = L-*F Ps (L y l v (K c -K),{ilL z ), 



(3) 



where F Ps is the scaling function, L is the linear system 
size, z the dynamical critical exponent, v the correlation 
length exponent, and (K c — K) cx (V c — V)/t is the dis- 
tance to the critical point. Thus if we plot p s L z as a 



function of V/t at fixed aspect ratio (3/L z , the curves 
for different system sizes should intersect at the critical 
point. The inset of Fig. [5] shows such a plot for z = 1 and 
P/L = 3/(4f) with a clear crossing point at (V/t) c ~ 19.8. 
The data is thus consistent with a continuous ST — X* 
transition with the dynamical exponent z = 1. To ob- 
tain the correlation length exponent v, we plot p s L as a 
function of [(V/t) c - V/tfL 1 ^ for different syste m sizes. 
It follows from Eq. ([3]) that the curves for different sys- 
tem sizes should collapse onto a universal curve F Ps for a 
properly chosen (V/t) c and v. In Fig. we show such a 
data collapse for (V/t) c = 19.80(2) and v = 0.67(5). The 
error bars are estimated from the stability of the collapse 
towards varying the parameters. To summarize, we find 
a continuous ST — X* transition with exponents z = 1 
and v = 0.67(5). We next examine the nature of the 
insulator X* . 



B. Characterizing T*\ Topological degeneracy and 
absence of broken symmetries 

In order to test whether the insulating phase of this 
model is a conventional broken symmetry state, we have 
studied correlation functions in X* . We have looked 
for signatures of diagonal (density), bond or plaquette 
ordering by studying the equal-time density and bond 
structure factors to check for different ordering patterns. 
Even for system sizes as large as 48 x 48 kagome unit 
cells and temperatures as low as T/J T i ng w 0.2, where 
Jring = t /V> we nn d no evidence of any Bragg peaks, 
or any ordering tendency, in these correlators. This ap- 
pears to rule out the possibility that X* is a conventional 
lattice symmetry broken state. Additional evidence for a 
uniform insulating phase comes from the fact that if the 
insulator had broken lattice symmetries it would not be 
smoothly connected to the high temperature paramag- 
net but be necessarily be separated from it by a thermal 
phase transition. However, we find no apparent signs of 
any thermal phase transition upon heating up from T = 
towards the uncorrelated high T paramagnet. 

For a system of bosons, momentum counting 



argument: 



,19,20 



show that an insulating state at half- 



filling could either be a conventional state with broken 
lattice symmetries or must necessarily have topological 
order which means the ground state degeneracy depends 
on the topology of the system. Since we have ruled out, 
as best as we can, the possibility that the insulating phase 
breaks lattice symmetries, we next examine the insulat- 
ing phase for signs of topological order. 

On a lattice with periodic boundary conditions in both 
lattice directions, the subspace of configurations with 
7iq = 3 on every hexagon of the kagome lattice is iden- 
tical to the Hilbert space of a triangular lattice quan- 
tum dimer model, with three dimers touching each site 
(identifying the hardcore boson with a dimer). It is well 
known that such quantum dimer models have well defined 
topological sectors, distinguished by the eigenvalues of a 
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FIG. 3: (color online). Average parities as a function of 
Monte Carlo time for two different system sizes and T = t/20 
in the insulator at V/t = 24. The start configuration has 
pi.2 = 0. The parities are finite because there is a small 
density of defects. The parities for L = 12 do not fluctuate 
indicating that the system does not change its topological 
sector whereas they strongly fluctuate for L — 1 indicating 
that the topological sector is changing. 



nonlocal operator, which do not mix under the dynam- 
ics generated by the Hamiltonian. In the context of our 
model, if we set t/V = and examine the classical ground 
states, which do respect uq — 3, these sectors corre- 
spond to having, for each lattice direction a\ t 2, an odd 
(or even) number of bosons on each row (or column) of 
the lattice (so called "parity sectors"). The row/column 
parities defined in this manner do not, however, specify 
topological sectors of model (p} since for any nonzero t, 
no matter how small, there will be a small density of de- 
fect hexagons^! with uq ^ 3 so that the row or column 
parity is not conserved under the Hamiltonian dynamics. 
How do we then look for topological sectors in the ground 
state of model (JTJ? 

We have checked in our QMC numerics, that if we 
start from a configuration which lies in the dimer sub- 
space with a certain row/column parity, then the QMC 
algorithm generates a small density of particle-hole pair 
defects in equilibrium, so that the quantum ground state 
no longer lies in the "dimer subspace" . However, for 
a large enough linear system size L at a given value of 
V/t (e.g., L > 10 at V/t = 26), our simulations with 
the longest accessible MC steps do not lead to any non- 
local boson moves which wind around the lattice, see 
Fig. [3J Thus the winding number identically vanishes, 
and each configuration in the simulation which lies in 
the dimer subspace belongs to the same parity sector as 
the initial configuration. The full ground state accessed 
by the QMC is, in this sense, perturbatively connected 
to the initial parity sector. In other words, we can start 
with the equilibrium QMC ground state and erase nearby 
particle-hole defects in pairs and obtain a state which lies 
entirely in the starting topological sector. It is in this 
sense that we can identify the four topological sectors 
even for model {!]), and we can continue to label them 



simply by the row/column parity of that component of 
the ground state wavefunction which lies in the dimer 
subspace. The four ground state wavefunctions can be 
thus be written as 



IV'oo! 



= 1^0, 



= \<1. 
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where IV'oh) denote the components of the wavefunctions 
that lie in the dimer subspace defined by parities a and 
b and \ip' ab ) denote the components of the wavefunctions 
that do not lie in the dimer subspace but are connected 
to \ipab) by short range hops of the bosons. It is clear 
that these wavefunctions are distinct eigenstates of our 
local Hamiltonian and they are not connected by local 
moves. For a topologically ordered insulator, the ground 
state energy should be identical in each of the four topo- 
logical sectors (on a torus) leading to a ground state de- 
generacy of four. We have computed the energy of the 
four ground states by starting our simulations from con- 
figurations in the dimer subspace lying in four different 
parity sectors. We find that they are equal within sta- 
tistical errors, which is strong evidence for topological 
order—. 



C. Vison correlations 

A second signature of Zi fractionalization is that vi- 
sons, which are gapped Z2 vortices in the effective field 
theory description, should have exponentially decaying 
spatial correlations. The spatial vison- vison correlation 
function is the expectation value of a string operator in 
terms of the spins. For the ring-exchange model (valid 
for V/t — > 00), it takes the forrn^: 



(4) 



k—i 



where |0) denotes the ground state and the product is 
along some path on the kagome lattice that contains an 
even number of sites, starts at site i, and ends at site j, 
making only "±60°" turns to the left or right, = 0, 1 
is the number of bosons at site k. The absolute value of 
the product in Eq. ^ is path-independent in the dimer 
subspace, and it is expected to decay exponentially in 
the topologically ordered phase. In model Jl} at finite 
V/t, ground state no longer lies entirely in the dimer 
subspace, but will mix in configurations with particle- 
hole defects. However, in the same manner as we have 
used the dimer subspace component of the wavefunction 
to define topological sectors, we can similarly use that 
wavefunction component to also compute the vison cor- 
relator. We have found that C vv (rjj), computed by essen- 
tially dropping all configurations containing particle-hole 
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FIG. 4: (color online). Energy per site E and compressibility 
k versus temperature for L = 24 (V/t=20.5). The energy rises 
exponentially at very low T (see inset) with a wide intermedi- 
ate plateau from T/t ~ 0.5 — 3 (see text). The compressibility 
is zero (within error bar) at low T, rising only at T ~ 2.5i 
once gapped charge (spinon) excitations become relevant. 

defects, decays exponentially in I*, again signaling topo- 
logical order in I* . At V/t = 20.5, we estimate a "decay 
length", £ = 1.43(5), which is comparable to its value at 
the RK point of the ring-exchange models, £ ~ 1, and to 
that found by EDi&, £ w 1.7, in the ring-exchange model 
with t>4 = 0. This exponential decay shows us that signif- 
icant local particle rearrangements and fluctuations are 
possible even in the insulating phase — such fluctuations 
are generated by the effective ring exchange dynamics of 
the bosons. 



III. HEATING THE INSULATOR: 
COOPERATIVE PARAMAGNETIC REGIME 
AND THE VISON GAP 

To provide further evidence for gapped vison excita- 
tions, we display the temperature dependence of the sys- 
tem energy per site and compressibility 

"*((?")') 

in Fig. |U Upon heating up from the ground state, the 
energy exhibits a two-step increase, with a distinct inter- 
mediate plateau. At the lowest temperature, the energy 
increases exponentially from the spin-liquid ground state 
as seen from the inset of Fig. |U The energy then reaches 
a plateau at a temperature T ~ J r ing- This plateau cor- 
responds to a regime where the system dominantly ex- 
plores configurations with uq — 3 on each hexagon. In 
spin language, it corresponds to a "cooperative paramag- 
net" . Upon heating further, the energy increases beyond 
its plateau value when the temperature becomes compa- 
rable to the local charge gap set by V. We confirm this 
by noting that there is a sharp increase in k at this higher 
temperature (also shown in Fig. 2]). 



Heating up from T = 0, we therefore identify the lowest 
energy excitations out of the ground state as coming from 
vison-pair excitations of the spin liquid (since the charge 
gap is much larger). The temperature dependence of the 
energy thus gives us a rough idea of the single vison gap; 
for V/t = 20.5, we estimate E v /t ~ 0.35(15). 

Note that the classical model with t = has a large 
entropy at T = arising from a macroscopic number 
of degenerate classical ground states. When one turns 
on a nonzero t, this degeneracy is lifted as the ground 
state becomes a quantum spin liquid which supports vi- 
son excitations [Z<i vortices). From the point of view 
of the spin liquid ground state, we can therefore view 
the large entropy density of the classical dimer state as 
arising from the large entropy density of multiple vison 
excitations as the system is heated. This means that the 
energy curve should have a finite temperature plateau at 
the energy level with the largest entropy density where, 
very roughly, half of all allowed visons get excited which 
contributes to a large configurational entropy of visons. 
This happens at an energy E s w Eq /2 as measured from 
zero classical ground state energy, where Eq is the quan- 
tum ground state energy. This is consistent with our 
numerical data. 



IV. HEATING THE SUPERFLUID: 
SUPERFLUID-NORMAL PHASE TRANSITION 

For V < V c , the system is in a superfluid ground state. 
Heating this superfluid leads to transition from a super- 
fluid phase to a normal liquid phase at finite tempera- 
tures. This transition is usually of a Kosterlitz-Thouless 
(KT) typ o 23 ' 24 that is driven by vortex unbinding. In 
principle, this transition can also be first order when the 
vortex core energy is small enoug h 25 i 26 . We find from our 
numerics that the thermal superfluid normal transition is 
a KT transition at small enough values of V/t. This goes 
along with the conventional wisdom. However, the tran- 
sition becomes first order roughly at V/t > 13(1). This 
appears to be a novel example of a discontinuous ther- 
mal superfluid-normal transition in a microscopic two- 
dimensional quantum model. 

A. QMC results for the SF-normal thermal 
transition 

In Fig. 03 we show the superfluid density p s and the 
system energy E as a function of temperature at V/t = 4. 
Both quantities exhibit smooth behavior. The super- 
fluid density should be discontinuous at the KT tran- 
sition with the universal jump 

a ZTkt 
Ap s = , 

where Tkt is the KT critical temperature. However, this 
discontinuity is approached only logarithmically as the 
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FIG. 5: (color online). The superfluid density p 8 and the 
energy E as a function of temperature at V/t = 16. 



system size increases. The RG equations — predict that 
in the vicinity of the KT transition the superfluid density 
scales as2&2L22, 



2T 



Pa = 



1 



F[(T-T KT )ln 2 (L/L ) 
2In(L/Lo) 



(5) 



where F is the scaling functions such that F(x) « 1 + 
0{x) for small values of x, L is the linear system size, and 
L is some length scale. If one plots (irps/T — 2) ln(i/T ) 
as a function of (T — Tkt) ln 2 (L/Lo) then it follows from 
Eq. O that the curves for different system sizes should 
collapse onto a universal curve F for a properly chosen 
Tkt and Lq. In addition, (irps/T — 2)ln(£/Xo) should 
take a value of 1 at the KT transition point. 

In Fig. [SI we show such a data collapse for Tkt = 
2.486i and Lq = 2. Thus we can conclude that the tran- 
sition at V/t = 4 is a KT transition. 

The situation is strikingly different for values of V/t 
roughly larger than 13. In Fig. we show the superfluid 
density and the system energy as a function of temper- 
ature at V/t = 16. Both the superfluid density and the 
energy jump suggesting that the transition is first or- 
der. As can be seen form Fig. [5J the distribution of the 
kinetic energy, Ek = —t((b\bj + H.c.)), close to the tran- 
sition has a clearly visible double peaked structure even 
for very small systems sizes indicating that the transi- 
tion is strongly first order. The double peaked structure 
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FIG. 6: (color online). Data collapse of the superfluid density 
ps for Tkt = 2.486£ and Lo — 2. Horizontal and vertical 
dashes lines are guides for the eye. 




FIG. 7: (color online). The superfluid density p s and the 
energy E as a function of temperature at V/t = 16. 



becomes more pronounced, leading to two well-separated 
peaks, as the systems size increases. We also observe 
hysteresis effects by crossing the transition point upon 
heating or cooling the system (not shown). This is also 
indicates that the transition is first order. 

For 13 < V/t < 17.5, the normal liquid just above T c 
has an energy which is nearly temperature independent 
for a range of temperatures. In this sense, it is analogous 
to the finite temperature plateau regime shown in Fig. 0] 
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FIG. 8: (color online). Distribution of the kinetic energy 
close to the transition for different systems sizes and slightly 
different temperatures (/3t = 1.67 for L = 6, j3t = 1.745 for 
L = 8, and (it = 1.77 for L = 12) at V/t = 16. A double 
peak structure is visible even for very small systems sizes and 
it becomes more pronounced, leading to two well-separated 
peaks, as the system size increases indicating a strongly first 
order transition. 

which is obtained upon heating the insulator. This nor- 
mal liquid can be well approximated by a classical liquid 
with no = 3 on each kagome hexagon, a so-called "coop- 
erative paramagnet" in spin language. However, when we 
increase V/t in the range 17.5 < V/t < V c /t, the normal 
liquid just above T c has increasingly significant quantum 
dynamics (as it lies in the quantum critical regime of the 
quantum phase transition) and, consequently, a progres- 
sively lower entropy than the classical liquid with no = 3 
on each kagome hexagon. The first order thermal transi- 
tion correspondingly gets weaker for V — > V c and eventu- 
ally merges into the continuous quantum phase transition 
at V c . 



B. Analysis of the thermal superfluid- normal 
transition 

We next turn to an analysis of the thermal superfluid- 
normal transition in order to shed some light on the ob- 
servation that it becomes first order rather being a con- 
tinuous Kosterlitz-Thouless transition when V/t is close 
to its quantum critical value. One qualitative way to 
understand this physics is to appeal to an entropy mis- 
match argument. At low T in the superfluid, the only 
relevant excitation is the superfluid sound mode. This 
has a velocity c s which is nonsingular at the z = 1 quan- 
tum phase transition into the insulator. The low tem- 
perature specific heat and entropy thus scale as ~ T 2 /c s 
deep in the superfluid. At large V/t the finite T normal 
state, if it lies in the "cooperative paramagnet" regime, 
is highly constrained and can be described very crudely 
by all possible classical configurations for which no = 3 
on each hexagon of the kagome lattice. This classical 
picture where we ignore all quantum fluctuations then 



predicts a large constant entropy S c \ in the normal state. 
Extrapolating from the low temperature superfluid state 
to the normal state thus leads to a large entropy mis- 
match, !§ c i — T c 2 /c s , at the transition when T c is small. 
Clearly, once T c decreases below a certain value, this en- 
tropy mismatch cannot be satisfied by the configurational 
entropy of a dilute gas of vortices. The superfluid to nor- 
mal transition must, at this stage, involve producing a 
large number of vortices over a small temperature inter- 
val to make up for this entropy mismatch, implying a 
large vortex fugacity. Such a large vortex fugacity is well 
known to modify the KT transition into a first order tran- 
sition. Of course this argument does not hold very close 
to the quantum critical point where the normal state just 
above T c is not a simple classical liquid with an entropy 
S c \, but is instead a quantum critical liquid. 

In what follows, we will provide a slightly different view 
of the thermal transition. It is well known that a Z2 frac- 
tionalized insulating ground state can be obtained from 
a superfluid ground state by condensing double vortices 
instead of single vortices. If the superfluid is proximate 
to such an exotic insulating ground state, as our other re- 
sults show, then the thermal excitations of such a super- 
fluid must include low lying double vortex excitations in 
addition to single vortices. Our aim here will be to infer 
the presence of these low energy double vortices from the 
observed first order superfluid-normal transition. We will 
therefore attempt to study the thermal SF-normal tran- 
sition taking into account both single and double vortices 
in a sine Gordon model. 



C. Phenomenological model 

Our discussion will use the classical sine Gordon model 
to describe the thermal phase transition from the super- 
fluid to the normal phase. To obtain this, we begin with 
the classical XY model written in vortex language 

(3H V = 2?r 2 K^2n r n r ,G{r -r') + J2(3E c {n r ) (6) 

rr' r 

where n r is the vortex number, G(r — r') is the vortex 
interaction (which is logarithmic at large distances), K 
is the superfluid stiffness normalized by the temperature, 
and E c (n T ) is the local core energy of a vortex with vor- 
ticity n r . We expect that at a fixed temperature T, the 
normalized stiffness K will decrease as we increase V/t. 
The hope is that such a classical description is adequate 
to qualitatively capture the physics of the thermal transi- 
tion with quantum fluctuations being important in fixing 
the parameters of this classical vortex Hamiltonian. 

We can go to k-space where the vortex interaction 
takes the simple form G(k) = 1/(4— 2cosfc;E — 2cosfc a ) for 
a 2D square lattice. Actually, the only important thing is 
that G(k — > 0) = 1/k 2 . The detailed form and the lattice 
geometry are unimportant. We can then do a Hubbard- 
Stratonovitch decoupling of the vortex interaction term, 
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writing 

so that the partition function takes the form 



where the coefficients A 



"i,i(g) r> _ "1.1(9) ^ _ 
i-ng ' -° 167rg ' ^ 1 



(»r} 

x e -E r ^K) 



(9) 



(10) 



Let us assign core energies (3E C (0) — 0, (3E C (1) = e\, 
f3E c (2) = e-i and /3E c (n > 2) = 00 to simplify the situ- 
ation. In that case, doing the sum over n r at each site 
leads to 

Z ~ j D4> r e-^"'^-**' )2/{ ^ K) (11) 

X ]^[(1 + Ul cos</> r + v 2 cos20 r ). (12) 
r 

where v p = 2c~ Cp . For small v p , we can re-exponentiate 
the cosine terms to get 

D^e-E^-^ 2 /^ (13) 

X e X)r(' UlCOS ^ r+1,2COS2 ^ r ) (14) 

Relabelling l/(8n 2 K) — > g/2, and in the continuum, 

Z SG = J D(j) r C~J d2r ^ W) 2 -»ico S *- 1 .2Cos2^ ^ 

The usual XY model corresponds to setting e 2 — 00, or 
1*2 = 0. Let us refer to these sine-Gordon models as SG2 
(with i>2 ^ 0, containing double vortices in addition to 
single vortices) and SG\ (with V2 = 0, with just single 
vortices). We can now study the RG flow equations for 
SG2, looking at the combined flow of <?, v\, 1*2. In terms of 
dimensionful quantities, which will be of use later while 
trying to understand the numerical phase diagram, we 
have K = p s /T (p s being the superfluid stiffness), so 
that 5 = T/(An 2 p s ) and v n = 2 exp[~(3E c (n)]. 

D. RG equations 



128^1 1 an d C-2 = "s'C^f • Here, we have defined 



dxx r 



; J {x). 



(19) 



The lower limit of this integral defining a m ,n{9) is set to 
unity — this reflects the short distance cutoff from the 
lattice spacing (which we have set to ~ 1/A) in doing 
various spatial integrals in this calculation. 

Physically, the second order correction in the flow 
equation for v\ may be viewed as arising from a single an- 
tivortex combining with a double vortex to give a single 
vortex. Similarly the second order correction to V2 can 
be viewed as arising from a double vortex splitting into 
two single vortices. For positive 0:1,1 (which is the case 
as we find numerically), the effect of a nonzero v 2 2> 1*1 
is that it enhances 1*1 due to the second order coupling 
Eq. (|16[) . The main observation we make is that integrat- 
ing out soft double vortices renormalizes the single vortex 
fugacity to be larger. 

We can simplify the above RG flow equations by work- 
ing around g ~ I/871". Over a range of g near this value, 
we approximate 01.1(5) « 5.55, 03,1(5) ~ 6.5<?, and 
03,4(5) ~ l-5g. Thus, A w 0.45, B w 0.11, d « 0.005, 
C2 ~ 0.02. With these simplifications, we have numeri- 
cally studied the above RG flow equations. 



E. Flows for vi <C 1 and V2 <C 1*1 

For initial V2(0) — and i*i(0) C 1, we recover the 
Kosterlitz Thouless RG flow. Namely, there is a line of 
fixed points (1*1,1*2,5) = (0,0,5*) with < 5* < l/8ir. 
The termination of this line of fixed points is the KT 
transition point at which the stiffness (normalized by the 
transition temperature) , K* — l/(47r 2 <?*) = 2/V, is a uni- 
versal number. Beyond this point, g flows to strong cou- 
pling, so that the superfluid stiffness K — > 0, signalling 
the non-superfluid phase. For a nonzero 1*2(0) <C 1*1(0), 
the above picture remains unchanged, i.e., the Kosterlitz- 
Thouless transition is unaffected by a small fugacity of 
double- vortices. 



The RG begins by writing </> = </>> + <^>< , where </*> has 
Fourier modes with momenta Ae~ M < q < A, while the 
slow field <j) < has the other Fourier components. A is the 
upper momentum cutoff. We fix it so that irA 2 = Air 2 , 
or A = 2^/tt. Integrating over the fast fields, we arrive at 
the following RG flow equations 

f| = „ l(2 -_i_ ) + A„« (16) 

£-«0-*i (17) 

% = dvf + Cvl (18) 



F. Flows for «i C 1 and V2 ~ 1 S> «i 

For large 1*2(0), we find that 1*2 (I) very quickly renor- 
malizes to small values since it is strongly irrelevant in 
the superfluid phase. However, in the initial stages of 
the RG flow, it significantly affects the flow of V\{£). 
While 1*1 (£) tends to decrease due to the first order term 
(2 — 1/4-715)1*1, this decrease is partially offset by the 
positive contribution from the second order term which 
couples 1*1 and 1*2- For our calculated couplings, in the 
regime where 1*1 is eventually irrelevant, it flows to zero 
more slowly for nonzero A. A comparison of the flows of 
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FIG. 9: Numerically computed flows for Vi,vz,g shown for 
two different choices of V2 (0) . We have fixed a small vi (0) = 
0.2 and 8Tr 2 g(0) — 0.5. The top panels shows the flows for 
A = 0.45 in the flow equations which is the result of our 
perturbative calculation. The bottom panels show the flows 
for an artifically large A = 4.5 to indicate that it leads to an 
enhancement of vi. 



V\ and g for small and large i>2(0) is shown in the upper 
panel of Fig. [5] 

A more striking result is obtained by artifically increas- 
ing the coefficient A in the flow equation for v\ . We find 
that increasing this to values larger than that given by 
the above perturbative calculation leads to a dramatic 
rise of V\ in the initial stages of the RG flow, due to 
coupling with v-i. This is illustrated in the lower panel 
of Fig. [51 where we have chosen A — 4.5 instead of our 
calculated result A rj 0.45. This modification of the flow 
equation is in a purely phcnomcnological spirit. Such an 
enhancement might be possible once we include the ef- 
fect of quantum fluctuations or higher winding number 
vortices which we have ignored, but which do become 
important near the superfluid-insulator quantum phase 
transition; however, this is beyond the scope of our per- 
turbative analysis. 

Fig. [9] is the central qualitative result of our RG cal- 
culation. Namely, the interconversion between double 
vortices and single vortices together with the low core 
energy (a large bare fugacity) for double vortex excita- 
tions can lead to a significantly enhanced fugacity for 
single vortices. At the same time, the double vortices are 
themselves irrelevant at long length scales. 

We will next use this strongly enhanced single vortex 
fugacity, obtained at intermediate length scales in the 
RG, to argue for a first order superfluid-normal thermal 
transition. Let us begin from the SGi sine-Gordon the- 
ory with ^(0) ^> wi(0), and follow the RG flows until we 
reach a fixed length scale £ where v%(£) <C fi(£)- At this 
stage, we can drop «2 altogether and study SG\ with only 
v\ ^ 0. We have shown that at some intermediate scale, 
can become large. In order to accommodate this 
large ^i(£) within the SG\ theory, the SGi action must 
be tuned to have a large bare fugacity for single vortices 
in other words, if we ignore double vortices (which 
we have shown is reasonable), the large fi(£) must be 



viewed as arising from a large v\ (0) . 



G. Analysis of the usual sine-Gordon model 

From our above RG analysis, we conclude that it makes 
sense to capture the effect of double vortices by studying 
the effect of a large bare v\ in the phase diagram of SG\ . 
We appeal to a variational method to study this following 
Ref. [3l|. The variational treatment for the sine-Gordon 
model SGi replaces the action 



S SG = I d 2 r 
by the variational action 
S = I ' d 2 v 



|(V</)) 2 - Wl cos< 



|(V0) 2 + -m<t> 2 



(20) 



(21) 



with the understanding that when g gets large enough, 
v\ will tend to pin the field <p to integer values leading to 
a mass for (j) fluctuations. 

At leading order, the variational free energy 



/var ~ f + ((SsG ~ S ))o- 



(22) 



is minimized for a mass m which satisfies the self- 
consistency condition 



.9 
m 
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m + 4ng 



l/8irg 



(23) 

(24) 



It is easy to show that for vi < 0.5, the mass is zero 
for g < <7kt = 1/87T while it increases continuously for 
g > 1/8-7T. This is the regime in which the continuous KT 
transition obtains in this model. At larger Vi > 0.5, the 
mass jumps discontinuously from zero to a nonzero value 
at a transition point g c < 1/Stt. This phase diagram is 
qualitatively illustrated in Fig. [TUTb) . 



H. Scenario for the SF-normal thermal transition 

We now appeal to the above results to understand the 
thermal transition from the SF-normal thermal transi- 
tion. Far from the quantum phase transition between 
the superfluid and the Z2 fractionalized insulator, the 
thermal disordering of the superfluid proceeds via the 
usual KT transition. As the zero temperature phase ap- 
proaches the quantum critical point however, the fugac- 
ity of double vortices at finite nonzero T becomes much 
larger than the fugacity of single vortices since the SF- 
X* quantum phase transition is driven by double-vortex 
condensation. We model this situation in our classical 
Hamiltonian by setting E c (2) — * near the quantum 
phase transition, and keeping a nonzero and large E c (l). 
We expect this effective classical description to be ade- 
quate so long as we are not too close to the quantum crit- 
ical point. Fig. [TUT a) illustrates the expected qualitative 
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phase close to the quantum phase transition point. This 
first order transition must be accompanied by a nonuni- 
versal jump in the superfluid stiffness which is larger than 
that predicted by KT theory. This is consistent with nu- 
merical observations presented in the previous sections. 



FIG. 10: (a): Assumed behavior of vortex core energies as a 
function of V/t in the effective classical model of vortices used 
to describe the thermal superfluid-normal phase transition. 
The solid lines are a schematic of the transition curve T C (V) 
as obtained from the SSE simulations. Light solid line shows 
the region where this superfluid-normal transition is of the 
Kosterlitz Thouless type. Dark solid line indicates where the 
transition becomes first order. The arrow indicates a constant 
T trajectory going through the SF- normal transition, (b): 
The corresponding phase diagram in the usual sine-Gordon 
model. The trajectory in the left panel translates into this 
modified trajectory in the (g,vi) plane. 



behavior of the single and double vortex core energies, as 
well as T c , as a function of V/t. 

Let us imagine a trajectory at fixed temperature but 
along increasing V/t. In this case, the bare superfluid 
stiffness decreases with increasing V/t, so that g is in- 
creasing monotonically. Let us assume for simplicity that 
the bare single vortex core energy E c {\) does not change 
with V/t. However E c (2) rapidly drops near the quan- 
tum phase transition upon increasing V/t. From our ear- 
lier discussion the large V2 then leads to an effectively 
larger v\. Thus, constant temperature trajectories in the 
(V/t,T/t) plane are expected to translate into trajecto- 
ries depicted in Fig.flQTb) in the (g, vi) plane. We argue 
that this may be responsible for the observed first or- 
der thermal transition from the superfluid to the normal 



V. CONCLUSION 

To summarize, we have studied the zero and finite tem- 
perature phase diagram of a model of hard core bosons 
with local interactions which exhibits a topologically or- 
dered Z2 insulating phase at zero temperature. In mag- 
netic language, this is equivalent to finding a quantum 
spin liquid phase of a S — 1/2 quantum magnet. We have 
presented a number of numerical results, and some ana- 
lytical arguments, in support of this identification. We 
have also studied the finite temperature phase diagram 
and identified a "cooperative paramagnet" regime, and 
seen that the superfluid to "cooperative paramagnet" 
transition is a first order transition rather than a BKT 
transition. Further work is needed to see if there is any 
connection between spin liquids found in simple model 
Hamiltonians, such as the one studied here, and the ex- 
perimentally observed spin liquids in quantum magnets 
on kagome lattices^. 
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